import math
from scipy import constants as const
import numpy as np


#electron density m-3
ne = np.array([1.61e20, 5.0e19])

#electron temperature KeV
Te = np.array([21.93, 7.0])

#DT average mass
mass_DT = 2.5 * 1.67262158e-27

#plasma minor radius
a = np.array([1.78, 0.8]) * 0.6

#plasma major radius
R0 = np.array([6.6, 1.6])

confinement_time = np.array([1.78,1.0])

#auxiliary heating power
P_aux = np.array([60.6 * 1.0e6, 27.0 * 1.0e6])

volume_plasma = const.pi * a * a * 2.0 * const.pi * R0

ratio_impurity = np.array([3.0e-4, 3.0e-4])

#impurity charge number
z_impurity = np.array([60.0, 10.0])

#DT effect charge number
z_effect_DT = np.array([1.0, 1.0])





#DT fusion reaction cross section
#1 barn = 1.0e-24 cm2
#ref: nuclear fusion reaction, figure 1.3
cross_section_fusion = np.array([0.5 * 1.0e-24 / 1.0e4, 5.0e-3 * 1.0e-24 / 1.0e4])

velocity_DT = np.sqrt(Te * 1.0e3 * const.e / mass_DT)

#DT fusion reaction coefficient
#ref: nuclear fusion reaction, figure 1.5
#reactivity_fusion = 5.0e-16 * 1.0e6
#the data of figure 1.5 may be wrong
reactivity_fusion = cross_section_fusion * velocity_DT
print("reactivity_fusion             :    ", reactivity_fusion)



energy_alhpa = 3.5 * 1.0e6 * const.e


#DT bremsstrahlung radiation power
P_DT = 5.0e-37 * z_effect_DT * ne * ne * np.sqrt(Te) * volume_plasma

#impurity radiation power
P_imp = 5.0e-37 * z_impurity * z_impurity * ne * (ne * ratio_impurity) * np.sqrt(Te) * volume_plasma

#background electron energy density
energy_density_e = ne * 1.5 * const.e * Te * 1.0e3

#plasma transport loss power, electron energy = ion energy
P_trans = (2.0 * energy_density_e / confinement_time) * volume_plasma


#alhpa particle self-heating power
P_alpha = (0.5 * ne) * (0.5 * ne) * reactivity_fusion * energy_alhpa * volume_plasma

print("==========================================")
print("volume_plasma     :    ", volume_plasma)
print("P_aux             :    ", P_aux/1.0e6)
print("P_alhpa           :    ", P_alpha/1.0e6)
print("P_trans           :    ", P_trans/1.0e6)
print("P_DT              :    ", P_DT/1.0e6)
print("P_imp             :    ", P_imp/1.0e6)
print("==========================================")
print(" ")
